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ABSTRACT 


In 1960 R. E. Kalman introduced a least squares concept 
that gives optimal estimates for the state of some dynamic 
Systems. Inclu@iea Ys a brieft historicad intmeductson 
leading to his work, a summary of his work, and the applica- 
tion of the theory to the sonobuoy reference system used on 
the S3A aircraft. Also, a tutorial development of certain 
M@Ponvivies used in the filter is presented in appendixes 


Pend B. 
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L. LNERODUCT UGH 


Aircraft launched sonobuoys form an important and 
growing role in ASW; however, it is necessary to know the 
pesition of the sonobuoys to fully utilize the information 
they transmit. In the past. accurate positions were avail- 
able only if the aircraft flew over the top of the buoy at 
low altitudes. For a large buoy pattern such a method 
consumes time which could be used for other aspects of the 
mession. The Navy’s newest ASW airnemaft, the S3A, will 
use a unique technique called the Sonobuoy Reference System 
(SRS) to estimate the sonobuoy's location. The SRS uses an 
interferometry principle coupled with a modern digital com- 
puter to update buoy positions. To increase the accuracy 
©f the estimated position, the data from the SRS are proc- 
essed through a Kalman filter. 

P The Kalman filter is a recent advance in estimation 
theory that has been widely used in the aerospace industry. 
ime term, filter, was advanced by electrical engineems to 
describe an input-output system. The Kalman filter was 
@esigned as a linear filter. A linear filter is a process 
minctech converts an input to an output in such a way that the 
output resulting from two simultaneous inputs is the addition 
of the two corresponding outputs [1]. 

Historically, estimation theory has been concerned with 
estimating the values of system parameters based on some 


observation of the system. The first major estimation 





technique was developed independently by A. M. Legendre and 
K. F. Gatiss. Although Legendre published his resulsge wim 

1806 while Gauss published his results in 1809, Gauss is 
generally credited with inventing the method of ledee Squares 
because he derived the method from fundamental principles 

and claimed to have used the technique to determine the 
parameters describing the orbital motion of celestial bodies 
in 1795. Gauss! parameters were the most probable value of 
the unknown quantities based on the data available. His 
definition of most probable values was; 

The most probable value of the desired parameters will 
be that in which the sum of squares of the differences be- 
tween the actually observed and computed values multiplied 
by numbers that measure the degree of precision is a mini- 
mum [2]. 

Following Gauss' book there was a long period where least 
Squares was the only widely accepted method for estimation. 
In 1912 R. A. Fisher published the method of maximum 
likelihood and thereby established the present theoretical 

background for estimation theory [3]. Maximum likelihood 
eevamation is cast in probability theory which requires 
mirormatvion on the probability distribution describing the 
object being measured. For an important class of problems, 
there is a direct correspondence between the maximum like- 
lihood and the least squares methods. Assuming the object 
being observed can be described by a gavissian distribuvien 
the two theories give identical results for the respective 


estimators. 


Following Fisher's work there was considerable effort 
am the field of communications to eliminate woancom moms 
from the signal. A. N. Kolmogorov and N. Wiener, in 1941 
and 1942, independently developed a linear minimum mean- 
square theory that provided the best separation of signal 
and noise [4,5]. The Wiener-Kolmogorov theory is a least 
squares process which differs from Gauss' theory primarily 
in that it allows the signal to vary with time with known 
Statistical qualities and the filter is required to be 
finear. The main problem associated with their filter is 
the difficulty in calculating the required parameters espe- 
seaolly for a complex system. Their results formed the 


Dasis for the subsequent work of R. E. Kalman. 


II. KALMAN THEORY 


Kalman's first paper on his filter theory was published 
in 1960 [6]. A summary of his results follows in a discrete 
time version which is the form used by the SRS. Notation 
used in this paper and appendixes is; capital letters and 
small letters represent vectors or matrices and their 
weomooncnius, respectively, subscripts refer to théewtime 
period being considered, carets indicate values that have 
been extrapolated through time and bars indicate that the 
Palues have been smoothed or weighted by some measurement 
mnormation. The equations, describing the system being 


Beseryed, must be of the ferm 


et 4) meee” Bua: ) 


X, is a nxl vector that describes the state of the system 
@o time i. Each component of X. is a state variable; there- 


fore, X. is called the state vector. is a nxn tran- 


pelea 
sition matrix that determines how the system propagates 
irom time i-l to i. U, is a nxl vector that is a gaussian 
random process which represents the uncertainties of the 


Ptacve equations in describing the true State. The statis-— 


tical properties of Us are 


EOUEa= 0 Reval ee 3) 
ELU,U,° J = 10 ions ales |. (3) 
E(U,U,"] = Yu, for all i. (4) 





pystem observations are given by 


Zell ay (5) 


i i al 


Z, is a pxl (o< n) observation vector. Bach componenues 


Zs is an observation on the system at time i that is related 
wo the state variables through M. a pxXn measurement matrix. 
Kalman theory requires M, to be such that the observations 
ee@eelinear combinations of the state variables. Ws is a 

mae vector Chat 1s a gaussian random process which represenvs 


me INeCeruciIny les Ol  UMewMedolieingsprocessces au Gime ei 


The statistical properties of Ws are 


ECW, J = 0 for 2 lumen (6) 
E(W,W. 0] = 0 for i # j, (7) 
E(W,W, J z My, for all i. (8) 


In addition, U and W are assumed independent therefore, 


E(W,U,"] for all i,j. (9) 


J 


i 
) 


Kalman proposed to estimate the state X. from the last 


estimate X_ and the observation Z, in such a way as to 


il 
minimize the mean-square error in the estimate. That is, 


7 \e 
m to be minimized. The selution to this problem Nas eeecenc 


known as the Kalman filter. The optimal estimate of the 


state at time i is 





~ 


X, = X. + K, (2, == M,X, J Claly) 
where X, = 95 4-184 Ci) 
The gain matrix Ky Pewee wecn ied by 

K, = C.M.‘(M.c.m.* + ¥ : (13) 

i ae ab clea sein Ws 

— = Jt 
where C. = 95 4-194 4-194 4-2 + ae (14) 
and ©, .=C,.-K, .M, .C (15) 


i-1 i-1l i-l i-l i-l 


where C, is the nxn covariance matrix of the state variables. 
in words, the optimal estimate of x. LS) 15 oe optimal estimate 
at time i1~l propagated to time i and then modified by the 
weiehting matrix Ks times the difference between the actual 
Bbeervabion’ at time Deana the .estimaved Value Of Sto 6s-cu— 
Peeion Se time i. Equations (11) through (15) are the basic 
equations of the Kalman filter. Initial requirements for 
the filter are Xs Cy» and the variances of U and W. 
Although the Kalman filter was designed for linear 
Syscems it can be modified for use with a nonlinear system 
@ewan the case of the SRS. Typically, the state equations 
and the observations are nonlinear and both must be "linear- 
ized" by making linear approximations of the equations to 
Ose in the filter. Such a method is reterred to 42s tio 
tended Kalman filter. The extended Kalman filter equations 


are analogous to the linear model equations [8]. Since the 


SRo utiligesslinesr statescquations we ume seqtiavionc sarc 
extended to include only nonlinearsobservations. The stave 


equation is 


x (16) 


eee eg 


The gaussian process U,_, was omitted from eq. 56) te 
facilitate later development. The omission is equivalent 
to assuming that the actual system has been modeled error 


free. The observation is 


Z. = h(X, ) + W 


‘ ’ (17) 


merce («sis a monlinear function of the state variables. 
The estimate of the state at time i is 


X, = x Se Zp n(X,)] (18) 


The gain matrix is 
(DAC Bee ee eo | (19) 


where D, is a nxp matrix of the partial derivatives of h(-), 


with respect to the state variables, evaluated at X, = X, 
oso, 


C, = 6 @ T 


ae ieee tet oe 


an 


eae Doc (21) 


oO ea et a 


OF 


i-1l 


ey 





The appreximactons Vsed- ine phere bende disk i iamietel |) vier 
can lead to unsatisfactory solutions due to divergence. 


Divergence occurs when C the covariance matrix of the 


ie 
state variables, approaches zero. This results in the 
weighting factors becoming unjustifiably small causing the 
incoming data to be weighted too little [7]. Thus, the 
solution becomes dominated by the time propagation of the 
state equations. One source of divergence is system equa- 
tions that do not properly model the physical system. This 
situation can occur if the linear approximations implicit 
in eq. (19) are not adequate in describing the nonlinear 
System. These approximations are analogous to discarding 
all the second degree and higher terms in a Taylor's series 
expansion of the nonlinear equations. Several approaches 


mer shaping the state covariance matrix to avoid filter 


divergence are discussed in [7]. 


1EAR 


171. SONSEUOY REFERENCE SYSTEM 


The SRS was designed to periodically update the posi- 
tion of the sonobuoys based on information received over an 
macenna system. Figure shows the location of the anvennas. 
These antennas provide three possible observations on the 
sonobuoy. 

Plone nnd darrect ion cosine’. 

2. Transverse direction cosine. 

3. ueeenee tmommairceratt to buoy. 

Renee initormation is available in limited conditions. ‘fhe 
normal mode of operation is with azimuth information. The 
aZimuth measurements are calculated from the phase differ- 
ence of the sonobuoy radio signals measured between the ends 
of antenna pairs. The antennas, with the shortest spacings, 
meerucea tO give Che general relative bearing to the buoy. 
Once this has been determined the computer selects a pair 
Seeanvennas that are further apart and refines the bearing. 
This is a repetitive process which continues until the anten- 
na pair is separated as far as possible. Antennas are loca- 
ted on the top of the fuselage and in the vertical tail to 
allow buoy signals to be received while the aircraft is 
manked in a turn. 

During normal operation, the SRS sequentially checks 
each of the 16 VHF channels to determine ii any Ssienals ave 


being received. If no signals are being received on a 


ii 
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channel that has a buoy assigned to it, the buoy location 
and the uncertainty of its location is updated by propaga- 
Ging the last estimate. should the ainrcrait be withinere— 
ceiving range of two or more buoys on the same VH} channel, 
mo data are processed. When the SRS determines that a sig- 
mere can besprocessed forsiniormation, the direction cosine 
is calculated and a command is sent to the aircraft naviga- 
tion system to determine the orientation and position of 
the aircraft. This information is passed to the filter, 
processed and then displayed to the crewmen on a cathode 


imeay tube. 
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LV. KAEMAN PIbVER FOR Dim soho 


The mathematical formulation of the problem utilizes 
mo cOOrdinave systems, the local and the aircratt sysvem. 
ie local coordinate system is a three=-dimensional, right- 
handed system with the x-axis pointing eastward, the X-Y 
plane target to the surface of the earth and the origin 
arbitrarily faxed. The aircraft attitude coordinate system 
meomesed tO reference the ailrcrafit orientation in roll, 
mcceh and heading. The aircraft attitude can be transformed 
miawe local coordinate values through a rotation matrix. 

The state equations for the system are based on the 
local coordinate system. The model of the buoy motion used 


wy the SRS is 


Xe4q = Xy ~ x. aie (22) 
Vag) = ¥y + ¥4 dt (25) 
| el (24) 
ees Oy (25) 


where dt is the time lapse between observations. In matrix 


notation they can be expressed as 


Xe ay = oX, (26) 


abe, 





where 


Se ee (27) 
"a x. 
a8 
V5 
and 
1 0 dat Oo 
co il ® dt (28) 
o oO WL 
1 -C€ OF 


m@emsubscript on the transition matrix, ¢, has been sup-= 
pressed since it is not changing with time. Also, the buoy 
motion as indicated by eqs. (22) thru (25) lies in a plane 
Darallel to the X=-Y plane of the local coordinate system. 
This does not imply that the ctrvature of the earth is not 
taken into account only that buoy motion in the Z-direction 
is negligible. 

The observation, Z, would be a 3x1 vector assuming the 
three possible observations were available simultaneously. 
This would require excessive computer time to invert the 
mesulting 3x3 matrix in eq. (19) for calculating the gain 
matrix. The SRS uses a sequential updating method to avoid 
this problem. Whenever multiple observations are available, 
they are processed independently thereby reducing W and Z 
to scalars. The observation for the direction cosine input 
is 


Z, = cos 6, + W, (29) 
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where W, iS Jesumed to meet the requirements sot Che scauco man 
noise in eq. (5). The observation is related to the state 


variables through the geometry in figure 2. 





hagoure 2 


oe is the aircraft heading, Cee are the coordinates of 


the aircraft position and (x ) are the coordinates of 


b>?" b 
the sonobuoy. Cos 6 can be expressed as 








cos 6 = cos (6, ~ $) (30) 
and since 
Yvry: 
cos 9 = 2 , (31) 
Xp ox 
Saige cas = (329 
where R = [(x,-x,)° +- (y,-y,)° 1’? (33) 
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rt LOllows thav 


Nae. 
weed -lyea: ; [svi ate ae 
cos 6, 5 cos’ 6), + 3 sin 6,,- (34) 
a i 
Since R. is nonlinear in the state variables Xy and sie 


and considering 86 and Ya, 25 Known QUaMehitece Gace 


joa Xai? 


observation can be expressed as in eq. (17) 


he Ao VGA te) le 35) 


where Xs and Ys refer to the buoy coordinates at time i. 
The Kalman filter equations are given by equations 
(16) thru (21) with by = 


6 for all k. The state variables 


covariance matrix would be 


CoCr) Stet) ens) eC) 
CO = 2) nei) 2) “EP (36) 
ce) ey) ces )  <ei@s)) 
e(4)" eqn)  ¢(9) eo) 
where e(1) = Var [x] (37) 
e(5) = Var [y] (38) 
e(8) = Var [x] (39) 
c(10)= Var [¥] (40) 


and the remaining elements are the respective covariances, 


eae. 


eCoyme= ov lyon. (41) 
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However, since the modeled system has no uncertainty built 
into it some method must be used to prevent divergence. 

To insure that C would not become too small c(8) and c(10) 
are replaced by Var [x] + (dAdt)* and Var [Ly] + (aAdt)* 
prior to the propagation of the C matrix in time. The 
value of dA is parameterized at .002 ft/sec and it repre- 
sents an acceleration uncertainty [9]. Examination of 

eq. (34) reveals that the partial sGosveesveR with respect 


femme velocity state vectors will be zero. Accordingly, 


D= (nh, hy 0 0) (42) 


oe and hy are demimed in Appendix 2). 


The required initial values for aircraft launched 


buoys are: 


1. State Vector. Coordinates are the aircraft drop 


coordinates and the component velocities are 
zero unless additional information is known. 
e. Covariance. Off diagonal elements are zero. 
c(1) and c(5) are H’ where H is the altitude of 
the aircraft in nautical miles. ec(8) and c(10) 


are (5 knots)“. 


3. Measurement variance. This quantity depends upon 


which sector relative to the aircraft the buoy 
is located, W = Gly ease except when the buoy is 
within 35° of either side of the rear longitu- 


dinal line when W = (1/2.4)°. 


I, 





V. SUGDABILITY Or TH KALMAN FILTER SOR thnks 


There are two requirements for a filter used in real- 
meme data processing. First, the storage and computational 
demands on the computer must be attainable and second, the 
accuracy requirements of the problem must be met or exceeded. 
The computer demands of a Kalman filter are high when com- 
pared to other methods [10]; however, the SRS uses several 
Techniques to reduce this demand. These techniques are 
directly related to the accuracy of the filter. Although the 
meamanr filter 1s optimal for a linear system, there is no 
meairancee that the real situation is modeled correctly. 

The errors introduced in modeling or in simplifying calcu- 
ime loons may result in marginal or unsatisfactory filter 
performance. An important result of the Kalman filter for 
Mempucer programming is that all the information from 

Z 2, <i Z LS ecomreai mete LT oe 


i-l 1 


mmocore the old observations. In addition, the recursive 


eliminating the need 


nature of the equations make them readily adaptable for use 


in acomputer. The basic steps in programming the filter 


are: 
1. Initialize X, C, and W. 
Extrapolate C (eq. 20). 
3. Compute weighting coefficients (eq. 19). 
4, Obtain measurement and estimate the state vector 


feds.) and i?) . 


rae 





5. Update C to reflect the last estimate (eq. 21). 

6. Retiteneto vstep se. 

Kalman's assumptions, of a linear system with additive 
Baussian noise, are approximated by the SRS. The nonlinear 
observation has been discussed previously. The gaussian 
distributed bearing error assumption may be significantly 
violated for certain flight conditions. Whenever the line 
of sight between the aircraft and sonobuoy does not change 
Significantly as compared to the data rate, the bearing 
errors may become biased due to multipath effects. Multi- 
pathing is caused by the radio signal reflecting off the 
aircraft surfaces prior to reaching an antenna. This re- 
sults in an apparent new position for the sonobuoy. Such 
Eetors would not be gaussian distributed and would not be 
eompensated for by the filter. In an effort to reduce the 
biasing, a 1/7 scale model was used to empirically determine 
constants to apply to the observations. When the SRS deter- 
mines what sector the sonobuoy lies in relative to the 
eercraft, the multipath constant for that sector is added 
we the observation to remove some of the biasing. Whether 
the constants from the scale model will be sufficient or 
whether constants for each aircraft will need to be deter- 
mined will have to be decided after the SKS is Plight Tested 


in early 1973. 


el 





APPENDIX A 


DERIVATION OF GAIN MATRIX COMPONENTS 


The components of the gain matrix, 


Ky. are calculated 


Girectly from the principle of minimum mean square error. 


Scalar equations are used instead of matrix equations to 


facilitate gaining insight into the various relationships 


and to emphasize the approximations. 


Only one weighting 


constant is derived; however, the remaining constants can be 


obtained in a similar manner. 


ihe equation, for the smoothed estimate of the x-coordinate 


oe thesbuoy, is 


“a 
— 


Xi41 ~ “i412 7 


Ky is to be calculated 
EL (X54. 7 X44) 
is a minimum. 
Now 
442° “i+] 
let 
E =) Se ac 
Xs ah a 


Ky C254) - hss) ¥447)) (Al) 
so that 

ei (A2) 
Keg Fey ley y PhO ay Vga I~ Xp4y 
Kot xydt + ky ENO Vga ya 


-h(x,+x,dt,y,+y,dt)] a ie x, dt 
(A3) 


(A oe 


and ‘2 = x, — xX. 


ae 





substitution into eq. (A3) yields 


€ Be ey Cheat: kK, [h(x +x, dt,y,ty,dt)tw,, 


Eee i i 1 


ie | eee ere a 5 FF 
“h(x, +x,dt,y,t+y,dt)] ; (A6) 
mm, addition , 


h(x, +x,dt,y,ty,dt) = Ot 81" gta ee 
aL 


= ss ae ane eae ical OL lier Pik Pe (AT) 


Using Taylor series expansion about the eh smooth 


estimate, the expression for h(*,°) can be approximated by 


h(x,+xdt,y,+ydt) + h, (x, +x aE Nice Tih Pike, 


A af 
+ hy (7, +kdt F447 4at)(-e, -2, at) (A8) 
where 
by(oy7) = SBE) and ny (++) = BE (A9) , (A10) 


Since the arguments for h(*,:) do not change and for ease of 
manipulation, the arguments of h(-,:) are suppressed and 


equation (A8) becomes 


Wet pee ae Md) + h (=e =e at) (All) 
eS Yo Ys VG 


aS 





Substituting eq. (All) into (A6) and squaring both sides 


yields 
' 2 ; 
a* = — + 2° at? + ral | + 2c, ema .te2e k,l ] 
ee sl i Aa : 4 44 ce (nea 
RIDES Bek if ] (Al2) 
Xs ALall 
where 
[ ie Weay + (e te at) °n’ + (ce te Bees) oa 
a eas Yi YG y 
- 2w,,,(e te dt)h_ - 2w,,.(e +e dtjh 
eee Xs Ky x 741 ae Oa y 
+ 2(e +e dt)h. (e te dt)h . (A13) 
se | arene yy 
Taking the expected value and defining 
os Ele’ ] 08 ECE ] 
al i at als 
| (Al4) 
i . oO | 2 
Cay Nae: Ele, e, J 5 = Elw,,/ feu 
aE a al (A20) 
ELw +1- ELws44e,. J ElWa4iex, 0 
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with similar expression for the y-quantitwtes yaelds 


of = of tr 0° at? oft Ki, [o°+(0% ABe =. deta at“)ne 
i+1 i 4 4) ee atl “ 
+ (of 5g Ge de+os dt“)he + 2(0% tos dt+o° at 
yn) Wane 9, y sy Ys 
{ 
2 2 2 2 2 
+ or* dt )h_h J] + 20°: dt - 2k_,h_ (ao~ to%* dt) 
XY 5 x y XX li a X, XX 
2 °C A 2 
ee a lel he +ao_*° dt) —- 2k..h dtla_: tae dat) 
Toy XY XY Asta XX- Xs 
-~ 2k.,h dt(o% +0%* dt) (A21) 
ay, may xy i 
i al 
Solving eq. (A21) for the Kay which minimizes 0° using 
to: 
standard calculus techniques results in 
2 2 2 2 O OQ 2 O 2 
h [co +20°+ dtto. dt°~] + h [oe to*)Gtt+o* + vt aac eed 
ie _ XO Xs XX_ Xs Y OXY, XY, XY; XY 5 
ie B 
(A22) 
where 
2 Cees ] O 2 2 °C 2 
Pe= co +h (co. +20.: dtto* dt ) + 2n h (o to* dtto «= dt 
x Xs XX 5 Xs XY XY, XVy XVs 
Nae dt) + h A 9° nay ae dttas ee). (A23) 
XY 5 YG a eg 


a> 





Similar results are obtained if the following matrices are 


substituted into eq. (19) and expanded: 





Mm oO erm Wee co. oo: 0 0 
Hee xy x 
i 5 nN 5 SL 
m1 0 Mel |r 0% mG 
C. = Yu VG V5 (A24) 
i OO TaMOPN [Mostics os Le 
7 KY 5 Pen 
OPo OC sali Vos. oc. oF. Oo” 
1, NA eon 
D, = (hy, hy 0 0) (A25) 
VW, = 6° ; (A26) 
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APPENDIX B 


DERIVATION OF PARTIAL DERIVATIVES 


Pigure Bl illustrates the geometry for the derivation 


of the partial derivatives of the measurement, cos @4. 





DIC dsh Ie aL 


a is the vector from the aircraft to the sonobuoy. 


ee is the vector representing the antenna pair making the 


observation. 
Pe, = (x,-x,)4 a (y,-y, 5 (ee) 
PB, = (x, -x, )i + (y, -y, DJ (B2) 
t e — ae bs 


The corresponding unit vectors@ere 


X, ~X 
Cagis ele 0. her: a 








an 





P OMe, Vy =—V 
P, = —1_“2 i + sliees j (BY) 
1g 1w 
where 
D 5 1/2 
R= [ (x, -x,) a Ca ] Cay) 


and Ky 1s the Penegn of the antenna pair. 











Hence, calc 
es X, -X Yi-y k 
7 foo see a oe t 
Bl > 1 ie. = 5 - . (BO) 
Cm ue 
Kt 
Let 
X, -X ae | 
t G ie iy 
oo = | Od (B7), (B8) 
al k o k 
vy t 
then 
ek Vi-y 
S ba ba 
ae + ¢, —=5 . (B9) 


Pieererentiating with respect to the state variables yiedds 








a(cos 0) _ “1 r Sah piste es Oa: 7 dR (B10) 
OX R 2 aX 
b | R b 
SmrIce , oR Xp 7X 
"| : (B11) 
Cc X, -X 
ONCOSmO) 2 al eee emia (B12) 
dX, R Ro 
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Similarly, 


s(cos 6) _ £2 JpoV%a 
3¥+, R 2 





fin matrix notation ; 


aX t 
b pe! 


ay 


(B13) 


(B14) 
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